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Abstract 


With  the  objective  of  estimating  the  remaining  forest  cover  on  Leyte  Island,  land  cover  classifications  using 
Landsat  4  TM,  Landsat  5  TM  and  Landsat  7  ETM+  were  performed  for  the  years  1989,  1994,  1996  and 
200 1 .  In  order  to  reduce  cloud  cover  and  shadow  obstruction,  a  total  of  37  images  were  processed  and  then 
fused  together.  A  Support  Vector  Machine  (SVM)  was  chosen  as  classifier.  The  accuracy  for  the  classified 
images  was  above  80%  for  each  image.  A  total  of  about  66,000  hectares  were  detected  as  deforested  in  the 
time  frame  from  1989  to  2001. 


Introduction 
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During  the  United  Nations  Framework  Convention  on  Climate  Change  (UNFCCC)  Conference  of  the 
Parties  (COP)  13  in  Bali  on  December  2007,  the  international  community  has  called  upon  countries  to 
explore  the  concept  of  reducing  emissions  from  deforestation  and  forest  degradation  (REDD)  as  a  new 
mechanism  to  combine  forest  protection  with  objectives  of  climate  protection,  biodiversity  conservation 
and  improvement  of  local  livelihoods.  The  project  “Climate-relevant  Modernization  of  Forest  Policy  and 
Piloting  of  REDD  Measures  in  the  Philippines,”  funded  under  the  International  Climate  Initiative  of  the 
German  Federal  Ministry  for  the  Environment,  Nature  Conservation  and  Nuclear  Safety  (BMU)  supports 
the  country’s  efforts  toward  forest  and  climate  protection,  and  the  development  of  appropriate  policy  and 
instruments.  The  successor  project  “National  REDD+  System  Philippines”  is  supporting  the  scaling  up  of 
technical  solutions  toward  a  national  system  for  measuring,  reporting  and  verification  (MRV)  of  REDD+ 
results.  The  present  study  has  been  prepared  during  an  internship  with  the  Deutsche  Gesellschaft  fur 
Internationale  Zusammenarbeit  (GIZ)  GmbIT  at  the  Philippines  Department  of  Environment  and  Natural 
Resources  (DENR)  under  the  National  REDD+  System  Project  from  March  to  September  2013. 

In  discussion  with  DENR,  Leyte  Island,  situated  in  the  Eastern  Yisayas  region  of  the  Philippines,  had 
been  selected  in  2009  as  a  site  for  piloting  of  REDD+  measures  under  the  Philippine  National  REDD- 
Plus  Strategy  (PNRPS).  The  current  processing  of  multispectral  optical  data  and  radar  data  to  detect 
deforestation  patterns  in  the  pilot  site  (see  Estomata,  2013a,  2013b)  needs  to  be  integrated  into  a  consistent 
workflow  covering  multiple  reporting  periods  (GIZ,  2012).  This  will  inform  a  forest  and  land  use  change 
assessment  for  the  establishment  of  data  on  preparing  reference  and  reference  emission  levels  for  REDD+ 
implementation  at  a  sub-national  level,  i.e.  Leyte  Island. 

The  Landsat  classification  performed  in  this  study  is  part  of  this  effort  to  achieve  highest  possible  temporal 
coverage  of  satellite  images  for  the  island  and  quantify  land  cover  and  use  changes.  It  can  therefore  be  seen  in 
conjunction  with  the  forest  resources  assessment  (FRA)  on  Leyte  Island  (DfS  Deutsche  Forstservice  GmbFf, 
20 1 3)  for  elaboration  of  an  internationally  compliant  system  for  REDD  MRV  as  outlined  in  the  conceptual 
approach  to  REDD+  MRV  in  the  Philippines  (Seifert-Granzin,  2013). 

A  total  of  37  images  from  1989  to  2001  were  acquired  from  the  United  States  Geological  Survey  (USGS). 
The  image  preprocessing  consisted  of  the  following  major  steps:  calibration,  atmospheric  correction  and 
cloud/-shadow  masking,  radiometric  normalization,  terrain  correction  and  image  fusion.  Subsequently,  a 
Support  Vector  Machine  (SVM)  Classification  was  performed  to  delineate  forest  vegetation,  non-forest 
vegetation  and  urban  area/bare  soil.  Finally,  all  four  classification  results  were  compared  in  order  to  detect 
deforestation.  Figure  1  displays  a  visualization  of  the  pursued  workflow.  Further  description  on  the 
integrated  data  products  and  the  single  processing  steps  is  given  in  the  following  chapters. 
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Figure  1.  Visualization  of  the  image  processing  workflow 

Green  boxes  represent  single  image  preprocessing;  turquoise  boxes  represent  processing  steps  incorporating  several  images  combined  to 
reduce  cloud  coverage;  blue  box  processing  steps  incorporate  multiple  cloud-reduced  image  mosaics  from  different  acquisition  years  for 
change  detection. 


1.1  Study  Area 

The  island  of  Leyte  is  located  in  the  Eastern  Visayas  region  of  the  Philippines  with  center  coordinates  of 
10°47’38”  N  124°48’39”  E.  It  expands  approximately  170  km  from  north  to  south  and  100  km  west  to  east, 
covering  an  area  of  about  7,368  km2.  Predominant  land  use  is  coconut  tree  plantation  with  trees  present  as 
pure  plantations  but,  more  commonly,  as  single  trees  in  forested  areas  in  different  proportions.  The  island  is 
traversed  by  a  mountain  range  with  peaks  up  to  578  m  high. 


1.2  Data 


For  classification  and  change  detection,  a  total  of  37  images  of  sensors  Landsat-4  Thematic  Mapper  (TM), 
Landsat-5  TM  and  Landsat-7  Enhanced  Thematic  Mapper  (ETM+)  were  used,  which  are  freely  available 
via  the  Earth  Explorer  service  of  the  United  States  Geological  Service  (USGS).  See  Tab.  8  for  a  full  list  of 
all  images. 

These  images  were  integrated  to  produce  four  cloud-reduced  mosaics,  which  incorporate  the  time  steps 
1989-90,  1993-94,  1995-96  and  2000-01. 

For  atmospheric  correction  (see  Section  2.2  for  details),  additional  data  products  were  included  to  enhance 
modeling  of  spatial  heterogeneity  of  the  atmosphere.  A  brief  description  is  given  in  the  following. 

Ozone  concentration  measurements  of  the  Total  Ozone  Mapping  Spectrometer  (TOMS),  onboard  the 
satellite  Earth  Probe,  were  included.  This  sensor  estimates  the  optical  depth  (transparency)  of  the  atmosphere 
by  correlating  measured  radiance  of  different  ultraviolet  wavelengths.  Due  to  their  different  atmosphere 
interaction  characteristics,  ozone  absorption  in  particular,  a  lesser  portion  of  the  sun’s  energy  is  transmitted 
through  the  atmosphere  and  reflected  back  to  be  perceived  by  the  sensor.  This  is  then  compared  to  direct 
measurements  of  sun’s  radiance  (solar  irradiance).  A  90%  daily  coverage  was  achieved  with  a  geometrical 
resolution  of  1.25°  longitude  and  1.00°  latitude  (McPeters  et  al.  1998:3,13;  Masek  et  al.  2006:69).  Refer  to 
McPeters  et  al.  (1998)  for  comprehensive  TOMS  sensor  and  algorithm  descriptions. 

Also  integrated  were  the  surface  pressure  and  precipitable  water  (water  vapor),  both  of  which  are  NCEP/ 
NCAR  Global  Reanalysis  Products  that  are  synergistic  multi-data  products  offered  by  the  US  National 
Center  for  Environmental  Prediction  (NCEP)  and  National  Center  for  Atmospheric  Research  (NCAR). 
NCEP  is  part  of  the  National  Oceanic  and  Atmospheric  Administration  (NOAA)  National  Weather  Service 
(NWS)  (NCEP/NCAR  2001).  Within  this,  data  from  surface  measurements,  radiosonde,  aircraft  and 
satellite  are  assimilated  together  with  supporting  modeling  approaches  (Kistler  et  al.  2001). 

Also,  Digital  Elevation  Model  (DEM)  GTOPO30  was  included.  This  multi-source  dataset  with  a  geometric 
resolution  of  30  arc  seconds  ( 1  km)  was  produced  in  1 996  as  a  joint  effort  of  several  international  organizations 
coordinated  by  the  US  Geological  Survey’s  Center  for  Earth  Resources  Observation  and  Science  (USGS 
EROS).  See  USGS  EROS  (2012)  for  reference. 
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Image  Preprocessing 
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In  order  to  combine  several  images  acquired  on  different  dates  throughout  the  year,  great  care  needs  to  be 
taken  in  radiometrically  matching  the  different  acquisitions.  The  preprocessing  steps  include  calibration 
(conversion  of  image  digital  numbers  to  top-of- atmosphere  [TOA]  radiance),  atmospheric  correction 
(conversion  of  TOA  radiance  to  surface  reflectance  [SR]),  cloud  and  cloud  shadow  masking,  radiometric 
normalization,  terrain  correction  and,  finally,  image  stitching.  These  steps,  including  the  necessary  scientific 
terms,  will  briefly  be  described  in  the  following.  All  steps,  except  preliminary  image  mosaicking  and  clipping 
as  well  as  the  final  image  stitching,  have  been  performed  with  the  Landsat  Ecosystem  Disturbance  Adaptive 
Processing  System  (LEDAPS)  software  (Masek  et  al.  2013).  The  open  source  code  is  freely  available  via  the 
Oak  Ridge  National  Laboratory  (ORNL)  Distributed  Active  Archive  Center  for  biochemical  dynamics 
(DAAC),  which  is  part  of  the  National  Aeronautics  and  Space  Administration  (NASA)  Earth  Observing 
System  Data  and  Information  System  (EOSDIS)  (Masek  et  al.  2013). 

2.1  Calibration 

The  conversion  of  the  raw  image’s  unit  less  digital  numbers  into  at-sensor  radiance  (ASR)  in  W*m-2*sr- 1  *  pm- 1 
(Watts  per  square  meter  per  steradian  per  micrometer)  is  calculated  by  using  calibration  factors  included  in 
a  metadata  text  file,  which  is  delivered  together  with  the  actual  image  data  of  Landsat  (NASA  201 1:1 17). 
Further,  optical  ASR  is  corrected  for  solar  zenith  (incidence  angle  of  the  Sun’s  radiation  approaching  the 
atmosphere),  Sun-Earth  distance  (mainly  depending  on  the  day  of  year)  and  solar  irradiance  (solar  radiance 
approaching  the  Earth,  W*m-2)  to  acquire  TOA  reflectance  (Masek  et  al.  2006:69).  TOA  reflectance,  also 
referred  to  as  planetary  reflectance  of  albedo,  is  given  in  percent.  Band  6  thermal  radiance  is  converted 
to  effective  at-satellite  temperature  (thermal  energy  exiting  the  atmosphere  expressed  in  Kelvin;  thermal 
energy:  electromagnetic  radiation  with  wavelengths  of  3-1 5  pm)  using  the  calibration  factors  recommended 
in  the  Landsat  Science  Data  Users  Elandbook  (NASA  201 1:120). 

2.2  Atmospheric  Correction 

The  conversion  of  TOA  reflectance  to  surface  reflectance  consisted  of  the  following  steps:  extraction  of 
aerosol  optical  thickness  (AOT)  with  the  dark,  dense  vegetation  method  (DDV)  (Kaufman  et  al.  1997)  and 
then  supplying  it  to  a  Second  Simulation  of  the  Satellite  Signal  in  the  Solar  Spectrum  (6S)  radiative  transfer 
algorithm  together  with  the  TOMS  ozone,  GTopo30  DEM,  NCEP  atmospheric  pressure  and  water  vapor 
data.  6S  models  atmospheric  scattering  and  gaseous  absorption  and  corrects  the  optical  TOA  reflectance 
for  it  in  order  to  acquire  surface  reflectance.  Within  this  procedure,  surface  pressure  data  are  used  together 
with  the  GTopo30  DEM  to  adjust  Rayleigh  scattering  to  local  conditions  (Vermote  et  al.  1997,  Masek 
et  al.  2006).  Rayleigh  scattering  describes  the  scattering  of  electromagnetic  radiation  on  particles  smaller 
than  its  wavelength  (e.g.  red  sky  sunset).  For  LEDAPS,  an  altered  version  of  6S  is  implemented,  which  was 
introduced  by  Kotchenova  et  al.  (2006)  (Vermote  &  Saleous  2007). 


2.3  Radiometric  Normalization 


When  combining  several  images  to  reduce  overall  cloud  cover,  their  radiometric  properties  must  be  as 
similar  as  possible.  All  images  in  consideration  should  therefore  be  seasonally  as  close  together  as  possible  to 
reduce  phenological  variability,  which  strongly  alters  surface  reflectance  of  vegetation.  Yet,  in  several  cases, 
temporal  coverage  of  land  surface  reflectance  was  chosen  over  radiometric  similarity. 

Other  causes  for  spectral  variability  are  cloud  cover,  cloud  shadow  obstruction,  change  of  land  use/land 
cover,  phenological  variations  and  uncorrected  atmospheric  alterations.  Due  to  the  nonlinearity  of  these 
impacts,  radiometric  normalization  between  images  should  not  be  performed  with  a  linear  approach.  In 
order  to  detect  change  between  two  image  acquisitions,  a  regularized  iteratively  reweighted  multivariate 
change  detection  (IR-MAD)  (Nielsen  2007)  was  performed.  IR-MAD  is  insensitive  to  linear  change  like 
gain  and  offset  (histogram  stretch),  and  most  atmospheric  correction  schemes,  and  therefore,  it  can  reliably 
detect  true  nontrivial  change  in  multispectral  data.  Once  pixels  with  least  change  have  been  determined, 
they  can  be  used  for  radiometrically  normalizing  the  images.  This  is  described  in  detail  by  Canty  &  Nielsen 
(2008).  Both  routines  are  implemented  within  the  iMAD  software,  whose  open-source  Python  code  was 
obtained  via  the  personal  homepage  of  Morton  J.  Canty  (2012).  An  exemplary  visualization  of  the  effect  of 
radiometric  normalization  can  be  seen  in  Figure  2. 

As  master  reference  for  the  normalization,  the  image  acquisitions  from  May  12  1989,  May  18  1994,  May 
23  1996  and  March  26  2001  were  chosen  for  each  mosaic.  The  choice  was  based  on  percentage  of  cloud 
coverage  and  overall  image  quality  mainly  depending  on  unmasked  cloud  shadow,  haze  and  thin  cirrus 
clouds. 


Figure  2.  Effect  of  radiometric  normalization  on  the  quality  of  a  stitched  image 

Left  uncorrected,  right  corrected  (5-4-3  cotor  composite). 


2.4  Cloud  and  Shadow  Masking 


The  LEDAPS  software  incorporates  the  Automated  Cloud  Cover  Assessment  (ACCA)  algorithm  as 
described  by  Irish  et  al.  (2006)  as  well  as  an  internal  surface  reflectance  based  mask  (SRBM).  ACCA  is  the 
standard  algorithm  used  by  the  Landsat  Project  Science  Office  (LPSO)  of  NASA’s  Goddard  Space  Flight 
Center  (GSFC)  to  estimate  percentage  of  cloud  cover  in  Fandsat  image  to  be  recorded  in  the  US  archive. 
The  algorithm  uses  the  spectral  bands  2  through  6  to  apply  26  decisions/filters  to  delineate  cloud  cover. 
A  detailed  description  of  the  algorithm’s  processing  steps  is  given  in  the  Fandsat-7  Science  Data  Users 
Flandbook  (NASA  2011:132-143).  As  reported  by  Vermote  &  Justice  (2010:21),  ACCA  outperforms 
SRBM  in  terms  of  commission  but  is  more  prone  to  omission  error.  In  this  study,  ACCA  was  found  to 
drastically  overestimate  cloud  cover;  wherefore,  it  was  excluded  from  further  investigations. 

For  the  purpose  of  this  study,  the  output  quality  assessment  (QA,  integrating  cloud  cover  and  shadow 
obstruction)  of  FEDAPS  was  converted  to  a  binary  mask.  The  QA  mask  is  delivered  as  1 6-bit  image  layer 
and  contains  bit-packed  quality  flags.  Converting  these  decimal  values  to  binary  format  reveals  unique 
combinations  of  binary  quality  flags  as  listed  in  Table  1.  For  instance,  the  decimal  value  6  converts  to  binary 
value  110,  which  reveals  that,  read  from  right  to  left,  this  pixel  contains  valid  data  and  was  masked  as  cloud 
by  the  ACCA  algorithm  but  not  by  the  FEDAPS-internal  cloud  masking  (as  bit  8  is  0). 

The  masked  (cloud  and  shadow  obstructed)  pixels  were  eroded  to  eliminate  small  pixel  objects  (noise).  The 
remaining  larger  objects  were  dilated  (expanded)  to  restore  them  to  their  original  size  including  a  buffer 
zone  around  masked  areas.  See  Figure  3  for  a  visualization  of  this  procedure. 


Table  1.  Surface  reflectance  bit  packed  quality  assessment  (QA)  flags 
(from  documentation  of  LEDAPS  model  product  [Masek  et  al.  2013]) 


bit  description 

0 

unused 

1 

valid  data  (Dyes) 

2 

ACCA  cloud  bit  (Dcloudy) 

3 

unused 

4 

ACCA  snow  mask 

5 

land  mask  based  on  DEM  (1=land) 

6 

Dark  Dense  Vegetation  (DDV) 

7 

unused 

8 

internal  cloud  mask  (1=cloudy) 

9 

cloud  shadow 

10 

internal  snow  mask 

11 

land/water  based  on  spectral  test 

12 

adjacent  cloud 

13-15 

unused 

2.5  Terrain  Correction 


Sun  illumination  varies  across  a  scene  due  to  terrain-induced  differences  in  orientation  towards  the  sun,  as 
well  as  inter-image  differences  caused  by  seasonal  changes  in  sun  elevation  and  aspect.  This  is  most  apparent 
in  dark  shadowed  valleys  and  very  light  slopes  in  (or  close  to)  zenith  to  the  sun.  In  order  to  reduce  this  effect 
and  normalize  reflectance  values  originating  from  the  same  land  cover,  the  C-Correction,  as  proposed  by 
Teillet  et  al.  (1982),  was  implemented.  The  approach  was  further  extended  to  NDVI-thresholded  correction, 
whereat  correlations  between  surface  reflectance  values  and  sun  illumination  are  evaluated  separately  for 
several  NDVI-thresholded  classes.  The  choice  for  C-Correction  was  based  on  the  findings  of  Richter  et  al. 
(2009),  who,  in  a  comprehensive  review,  found  this  method  to  be  superior  to  other  prevalent  approaches. 


Figure  3.  Performance  of  cloud  masking 

From  left  to  right:  original  image  (4-5-7  composite),  LEDAPS  quality  assessment,  fusion  of  manipulated  mask  and  image.  Error  of  commission 
was  neglected  in  favor  of  low  error  of  omission. 


An  exemplary  visualization  of  the  algorithm’s  performance  can  be  seen  in  Figure  4.  Although  the  method 
performs  well  on  little  to  moderate  slopes,  shadow  in  terrain  with  steep  slopes  is  missed,  and  reflectance 
of  sun-zenith  aligned  slopes  is  further  enhanced.  This  is  primarily  related  to  the  lower  90-m  geometrical 
resolution  of  the  SRTM  DEM  and  is  at  this  moment  inevitable.  These  errors  are  partly  compensated  for  by 
a  later  performed  Minimum  Noise  Fraction  rotation  (MNF). 


Figure  4.  Image  subset  (4-5-7)  showing  the  performance  of  the  C-Correction  for  solar  illumination 

White  pixels  were  masked  for  cloud  and  shadow  obstruction. 
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Image  Fusion 
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The  master  images  used  to  start  the  fusion  were  the  same  images  used  for  the  radiometric  normalization  in 
Section  2.3.  The  stitching  was  then  performed  in  a  straight  forward  fashion,  where  if  a  pixel  of  the  master  or 
preliminary  step  fusion  image  is  masked  as  cloud  or  cloud  shadow  and  the  corresponding  pixel  of  the  slave 
image  is  clear,  this  clear  pixel  is  used  to  synthesize  the  new  image.  See  Figure  5  for  an  exemplary  visualization 
of  the  performance. 
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Figure  5.  Subset  of  Eastern  Leyte  showing  the  performance  of  the  cloud  stitching,  color  composite 
4-5-7  (NIR-SWIR1-SWIR2) 

From  left  to  right:  master  image  acquired  on  Mar.  26,  2001,  stepwise  combination  with  image  of  date  Dec.  4,  2000;  Feb.  22,  2001;  Apr.  11, 
2001.  Ctouds  and  their  shadows  are  masked  white,  whereas  water  is  black. 


Classification  and  Validation 


As  classifying  algorithm,  a  Support  Vector  Machine  (SVM)  was  chosen  as  implemented  in  the  ENVI  5.0 
software  (refer  to  CHANG  &  LIN,  2001  and  HSU  et  ah,  2007  for  algorithm  descriptions).  This  choice  was 
based  on  comparisons  with  other  classifying  algorithms.  SVM  showed  to  be  least  prone  to  errors  originating 
from  spectral  inconsistency  caused  by  insufficient  radiometric  normalization  and  cloud  shadow  masking. 
Training  subsets  were  chosen  across  the  island  and  validated  using  high-resolution  Google  Earth  imagery. 

For  enhancing  separability  for  choosing  training  subsets  and  the  classifier,  a  Minimum  Noise  Fraction 
(MNF)  Forward  Transformation  was  performed  for  Landsat  bands  1,  2,  3,  4,  5  and  7.  The  six  resulting 
layers  were  then,  after  choosing  representative  samples  of  the  required  classes,  supplied  to  the  SVM  routine. 
The  MNF  transformation  separates  between-band  spectral  correlation  in  multispectral  data  according  to 
the  Principal  Component  Analysis  (PCA)  principle  and  orders  the  resulting  components  by  their  signal-to- 
noise  ratio  (SNR)  (Green  et  al.  1988:65ffi).  Concerning  the  selection  of  training  samples,  this  procedure 
was  particularly  beneficial  as  inter-image  radiometric  inconsistency  was  enhanced  and  could  be  sampled 
appropriately.  Mentioning  the  insufficient  illumination  correction  in  mountainous  areas,  in  some  cases,  this 
effect  was  found  to  be  partly  eliminated  in  the  noisiest  MNF  components.  Yet,  this  effect  still  introduced 
the  largest  error  for  the  forest  class. 

Figure  6  displays  the  spectral  variability  in  the  used  image  channels  for  the  most  critical  classes  (forest, 
transition  and  palm)  as  boxplots. 

For  assessing  the  accuracy  or  quality  of  the  classification,  another  set  of  test  areas  was  chosen  from  ground 
truth  acquired  in  July  2013  and  comparison  with  Google  Earth,  and  then  checked  against  the  classification 
of  the  area.  The  test  samples  were  chosen  in  the  Landsat  image,  as  direct  definition  from  ground  truth 
polygons  and  the  Google  Earth  image  would  introduce  additional  error  due  to  the  large  difference  of  the 
acquisition  dates.  The  results  for  each  year  are  displayed  in  Tables  2  to  5. 

Overall  accuracy  was  found  to  be  between  81  and  90%.  The  Kappa  coefficients  were  0.74,  0.87,  0.83  and 
0.77  for  the  years  1989,  1993,  1995  and  2000,  respectively.  Largest  errors  are  encountered  by  delineating 
forest.  The  Urban/Bare  class  was  matched  very  well  due  to  its  strong  contrast  to  the  vegetated  areas. 
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Figure  6.  Spectral  variability  of  forest,  palm  and  transitional  zone  (mixed  pixels)  reflectance  and  MNF 
eigenvalues  for  selected  subsets 

The  red  line  represents  the  mean  value;  the  blue  box  delineates  variance  and  the  whiskers  mark  minimum  and  maximum  values. 


Table  2.  1989-1990  accuracy  assessment 


Class 

Forest 

Palm 

Mix 

Agri/Bare 

Total 

Commission 

Forest 

348 

20 

3 

1 

372 

93.55 

Palm 

53 

565 

87 

0 

705 

80.14 

Mix 

164 

63 

551 

1 

779 

70.73 

Agri 

91 

8 

20 

664 

783 

84.8 

Total 

656 

656 

661 

666 

2639 

Omission 

53.05 

86.13 

83.36 

99.7 

80.64 

Top-down:  test  samples  vs.  classification  (omission  or  producer's  accuracy);  left-right:  classification 
vs.  test  samples  (commission  or  user’s  accuracy). 


Table  3.  1993-1994  accuracy  assessment 


Class 

Forest 

Palm 

Mix 

Agri/Bare 

Total 

Commission 

Forest 

580 

16 

28 

0 

624 

92.95 

Palm 

66 

602 

18 

2 

688 

87.5 

Mix 

51 

74 

610 

2 

737 

82.77 

Agri 

0 

1 

2 

647 

650 

99.54 

Total 

697 

693 

658 

651 

2699 

Omission 

83.21 

86.87 

92.71 

99.39 

90.37 

Table  4.  1995-1996  accuracy 

assessment 

Class 

Forest 

Palm 

Mix 

Agri/Bare 

Total 

Commission 

Forest 

528 

38 

44 

0 

610 

86.56 

Palm 

134 

631 

22 

1 

788 

80.08 

Mix 

35 

13 

555 

0 

603 

92.04 

Agri 

8 

3 

40 

666 

717 

92.89 

Total 

705 

685 

661 

667 

2718 

Omission 

74.89 

92.12 

83.96 

99.85 

87.56 

Table  5.  2000-2001  accuracy 

assessment 

Class 

Forest 

Palm 

Mix 

Agri/Bare 

Total 

Commission 

Forest 

581 

65 

72 

6 

724 

80.25 

Palm 

41 

541 

10 

3 

595 

90.92 

Mix 

108 

130 

553 

1 

792 

69.82 

Agri 

11 

8 

54 

737 

810 

90.99 

Total 

741 

744 

689 

747 

2921 

Omission 

78.41 

72.72 

80.26 

98.66 

82.57 
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5.1  Preliminary  Agriculture  Mask 

For  detecting  true  change  in  land  cover/ use,  areas  of  agricultural  use  need  to  be  identified.  Due  to  their 
strong  phenological  variation  and  change  to  bare  soil  across  seasons,  areas  of  true  change  can  easily  be 
overestimated.  With  the  aim  of  reducing  this  error,  the  NDVI  was  derived  for  each  image  acquisition  of  one 
time  frame  (year),  maximum,  minimum  and  their  difference  calculated  and  thresholded  to  identify  areas 
that  have  undergone  large  changes  in  NDVI  following  conditions. 


NDVIdiff  >  0.1;  NDVImin  <  0.6, 


where  NDVId.^  is  the  difference  of  maximum  and  minimum  NDVI  value  for  each  pixel,  along  the  specified 
time  frame,  and  NDVImjn  is  the  minimum  value.  Thresholding  was  particularly  necessary  for  excluding 
larger  forest  areas  that  expressed  large  NDVI  changes  across  images  due  to  insufficient  correction  or  masking 
of  atmospheric  effects.  See  Figure  7  for  a  visualization  of  the  method’s  performance. 


Figure  7.  Masking  of  agricultural  area 

From  left  to  right:  5-4-3  composite  image  subset,  differential  NDVI,  thresholded  mask. 


5.2  Post-Classification  Change  Detection 


Change  is  detected  in  the  manner  of  cross-comparing  all  classification  results  of  the  time  frame.  If  one  pixel 
is  classified  as  forest  or  mix  of  natural  trees  and  palm  trees,  but  is  classified  as  agriculture  or  bare  soil  in 
the  following  time  step,  it  is  marked  as  deforested.  This  condition  then  only  holds  true  if  the  pixel  was  not 
classified  as  bare  soil  or  agriculture  in  any  time  step  before  the  forest  classification.  Lastly,  all  pixels  identified 
as  deforested  are  cross-checked  with  the  preliminary  agriculture  mask  discussed  in  Section  5.1  and  are 
excluded  if  they  have  been  identified  as  agriculture. 


Table  6.  Area  statistics  of  classification  results  in  hectare 


[ha] 

1989-1990 

1993-1994 

1995-1996 

2000-2001 

Fusion 

Forest 

50713.02 

106204.32 

111518.46 

1  10662.65 

146829.15 

Mix 

106168.05 

100042.38 

123214.41 

134953.65 

139180.68 

Palm 

50110.38 

70576.65 

122552.73 

90574.47 

87022.98 

Agriculture 

1  09126.35 

74288.43 

190131.12 

60366.42 

213587.55 

Bare/Urban 

217546.47 

77249.70 

53858.34 

156424.23 

53850.87 

Masked 

188020.08 

274581.27 

99653.85 

141548.40 

15324.30 

Defol 

0.00 

5168.34 

10436.13 

13091.49 

28502.37 

Defo2 

0.00 

13668.57 

10422.72 

14166.45 

37484.01 

Defo  1  and  Defo2  are  Deforestation  of  Forest  and  Mix,  respectively. 


Table  7.  Area  statistics  of  classification  results  in  percent 


[%] 

1989-1990 

1993-1994 

1995-1996 

2000-2001 

Fusion 

Forest 

7.03 

14.71 

15.45 

15.33 

20.34 

Mix 

14.71 

13.86 

17.07 

18.70 

19.28 

Palm 

6.94 

9.78 

16.98 

12.55 

12.06 

Agriculture 

15.12 

10.29 

26.34 

8.36 

29.59 

Bare/Urban 

30.14 

10.70 

7.46 

21.67 

7.46 

Masked 

26.05 

38.04 

13.81 

19.61 

2.12 

Defol 

0.00 

0.72 

1.45 

1.81 

3.95 

Defo2 

0.00 

1.89 

1.44 

1.96 

5.19 
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Forest  cover  and  palm  land  use  were  reliably  well  classified  within  the  scope  of  this  study.  Yet,  within  the 
processing  chain,  several  steps  can  be  improved  to  further  increase  the  quality  of  the  final  product.  Some 
of  the  flaws  will  be  described  in  the  following  to  outline  potential  for  further  investigations.  For  this,  not 
the  theoretical  basis  of  particular  algorithms  implemented  elsewhere  is  questioned  nor  validated,  as  this 
was  out  of  scope  of  this  study.  Rather,  focus  is  put  on  reliability  of  the  processing  chain  and  applicability  to 
succeeding  studies. 
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Figure  8.  Land  cover/use  and  deforestation  map  of  Leyte  Island  for  the  years  1989  to  2001 


6.1  Cloud  and  Cloud  Shadow  Masking 


The  algorithm  was  found  to  overestimate  cloud  shadow  in  various  situations,  particularly  mistaking  dark 
bare  soil  and  riverbeds  for  shadow.  Masked  clouds  were  not  precisely  outlined  in  many  cases  (over-  and 
underestimation),  which  demanded  further  editing  of  the  mask  through  eroding  and  dilating.  Cloud  and 
cloud  shadow  obstruction  of  the  final  product  could  be  decreased  by  implementing  a  more  precise  and 
therefore  less  generous  masking,  declaring  a  larger  area  as  free  of  clouds  and  shadow.  An  implementation 
of  a  refined  procedure  is  considered  rather  time  consuming,  as  the  applicability  to  other  atmospheric  and 
seasonal  situations  needs  to  be  comprehensively  evaluated.  The  LEDAPS  QA  mask  was  found  to  omit 
several  bit  combinations.  Bits  9  (cloud  shadow)  and  12  (adjacent  cloud)  never  exist  together,  which 
complicates  manipulations  on  singular  bitmasks.  An  image  object  segmentation  approach  is  suggested  for 
further  refining.  A  more  critical  difficulty  presents  the  insufficient  masking  of  occasionally  encountered 
large  cloud  shadows.  This  case  is  most  noticeable  in  areas  of  mixed  trees  or  palm  plantations,  which  are 
accordingly  mistaken  for  forest. 


6.2  Radiometric  Normalization 

Although  the  MAD  procedure  is  significantly  enhancing  the  quality  of  the  final  product,  visual  inspection 
revealed  room  for  improvement.  Other  methods  are  suggested  in  literature,  which  have  been  shown  to 
perform  better.  Refer  to  Helmer  &  Ruefenacht  (2007)  for  a  comparative  study  on  different  approaches. 
Here,  a  regression  tree  method  is  suggested,  which  can  be  implemented  in  open  source  versions.  Yet, 
documentation  on  the  application  of  this  code  (written  in  R)  is  poorly  documented  and  requires  further 
research  beyond  this  study.  However,  by  carefully  selecting  training  subsets  on  the  MNF  transformed  image, 
poorly  normalized  image  mosaic  patches  can  be  accounted  for.  Extreme  cases  of  insufficient  normalization 
were  only  found  with  images  of  extreme  cloud  obstruction,  incapacitating  the  MAD  routine  of  finding  an 
adequate  set  of  unchanged  samples  for  radiometric  correlation  computations.  This  case  could  be  accounted 
for  by  redesigning  the  MAD  and  radiometric  normalization  algorithms  to  more  efficiently  use  available 
memory.  Currently,  the  whole  image  is  read  at  once,  therefore  only  allowing  smaller  subset  processing. 
Implementing  band-wise  array,  computing  in  combination  with  water  masking,  would  enable  sample 
selection  across  the  whole  island  for  specifying  the  linear  regression. 

Further  research  is  also  suggested  to  determine  the  effect  of  image  acquisition  date  and  therefore  the 
phenological  variability  on  the  quality  of  normalization. 


6.3  Classification  and  Validation 

The  SYM  method  has  been  proven  to  be  very  robust  and  flexible  compared  to  other  approaches  and  is 
gaining  wide  acceptance  in  the  field  of  research,  see  for  e.g.  Townshend  et  al.  (2012:13),  who  used  this 
method  to  create  global  forest  maps  with  Landsat  data.  In  this  study,  SVM  was  also  observed  to  be  less 
mistaken  by  spectral  changes  related  to  unmasked  shadow. 

The  largest  uncertainty  when  classifying  images  lies  in  the  selection  of  training  samples.  Reference  satellite 
data,  with  beneficial  resolution,  spectral,  temporal  and  radiometric,  might  not  be  available.  Ground  truth 
campaigns  are  time  consuming  but  necessary  to  gain  comprehensive  understanding  of  land  cover.  Yet,  cross¬ 
comparison  with  images  acquired  in  different  years  is  to  be  performed  with  caution,  as  essential  change 
might  occur  in  small  time  steps. 

Separability  of  forest,  mix  and  palm  was  observed  to  strongly  depend  on  SWIR  reflectance,  which  is  closely 
related  to  moisture.  The  effect  of  seasonal  moisture  changes  on  separability  of  these  classes  and  is  yet  to  be 
assessed  in  a  wider  scope.  This  might,  if  changes  are  little,  allow  for  selecting  images  across  the  whole  year 
for  mosaicking. 


6.4  Change  Detection 


As  inherent  to  post-classification  change  detection,  the  accuracy  of  the  result  is  strongly  dependent  on  the 
accuracy  of  the  single  classification  results  and  never  better  than  those.  Yet,  in  the  scope  of  this  study,  it  was 
assumed  to  offer  the  best  compromise.  Methods  like  Change  Vector  Analysis  (CVA,  Malila  1980)  or  MAD 
(Canty  &  Nielsen  2008)  are  well  suited  for  detecting  spectral  change,  yet  the  nature  of  change  remains 
concealed.  Moreover,  change  can  have  numerous  causes,  with  unmasked  cloud  shadow  and  changing 
illumination  in  steep  mountainous  areas  being  the  most  conceiving.  In  this  study,  change  was  therefore 
limited  to  the  most  distinctive  case,  where  a  pixel  is  marked  as  deforested  if  it  was  forest  or  mix,  and  then  is 
bare  soil  or  agriculture.  The  spectral  signatures  of  these  classes  are  very  well  separable  and  unlikely  mistaken 
by  the  classifier,  even  though  the  actual  bare  soil  area  might  be  much  smaller  than  the  pixel.  The  largest  error 
in  this  procedure  is  introduced  by  the  ambiguity  of  the  agriculture  class,  which  frequently  undergoes  these 
changes.  The  accuracy  of  identified  change  is  therefore  strongly  dependent  on  the  temporal  resolution  of 
the  pixel.  More  samples  of  each  pixel  are  available  per  year  and  make  the  identified  change  less  ambiguous. 

Accordingly,  it  is  necessary  to  observe  that  the  accuracy  of  the  identified  deforestation  varies  across  the 
image,  depending  on  the  sampling  of  each  pixel.  This  is  dependent  on  the  number  of  cloud-free  acquisitions 
per  time  step  and  the  time  step  itself.  Later  observed  change  is  more  accurate,  as  more  cross-checks  with 
earlier  images  can  be  incorporated. 
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The  pursued  procedure  of  integrating  several  Landsat  images  to  classify  mosaicked  images  and  then  detecting 
deforestation  changes  from  1989  to  2001  was  found  to  be  very  reliable  in  encountering  difficulties  as  cloud 
cover  and  spectral  variability  gave  very  satisfactory  robust  classification  results.  Calibration,  atmospheric 
correction,  radiometric  normalization  and  cloud  masking  enabled  the  creation  of  homogenous  image 
mosaics  fused  together  of  images  acquired  on  various  dates  and  conditions.  Throughout  this  paper,  several 
indications  about  further  improving  the  workflow  are  given.  The  most  striking  ones  are  radiometric 
normalization  and  a  more  complex  inter-dependency  change  detection  approach. 

This  approach  will  further  improve  if  more  images  are  incorporated,  particularly  in  the  case  of  change 
detection,  as  this  is  specified  by  cross-checking  with  different  image  acquisitions.  The  biggest  limiting 
factor  in  this  analysis  is  cloud  and  cloud  shadow  obstruction,  which  strongly  limit  the  study  area.  The 
established  workflow  enables  the  analyst  to  easily  expand  the  scope  of  analyzed  image  frames  through  the 
use  of  primarily  automatized  procedures. 

With  the  perspective  of  installing  a  nationwide  satellite  MRV  system,  this  study  is  seen  as  an  important 
framework  for  tracking  medium  geometric  resolution  land  cover  and  land  use  changes,  which  can  easily 
be  extended  to  different  areas  and  time  frames  as  long  as  a  sufficient  number  of  images  is  available.  This  is 
not  always  ensured  due  to  very  high  percentage  of  cloud  obstruction  and  Landsat  technical  problems  like 
inconsistent  spatial  coverage  for  Landsat  7.  Use  of  different  datasets,  like  Synthetic  Aperture  Radar  (SAR) 
and  higher  geometric  resolution  optical  images,  is  strongly  recommended  to  enhance  temporal  accuracy  and 
small-scale  detail  of  deforestation  mapping. 


References 


Canty,  M.J.  (2012):  Personal  Homepage.  <mcanty.homepage.t-online.de>  (Status:  2012-08-08)  (Access: 
2013-05-15). 

Chang,  C.-C.  &  C.-J-  Lin  (2001):  LIBSVM:  a  Library  for  Support  Vector  Machines.  National  Taiwan 
University,  Department  of  Computer  Science  &  Information  Engineering  NTU  CSIE  online. 
<http://www.csie.ntu.edu.tw/-cjlin/papers/csvm.  pdf>. 

DfS  Deutsche  Forstservice  GmbH  (2013):  “Forest  Carbon  Baseline  Study  Leyte.  Final  draft  report.  GIZ 
German  Technical  Cooperation  with  the  Philippines.  Climate  Relevant  Forest  Policy  and  Piloting 
of  REDD.” 

Estomata,  M.T.  (2013A):  “ALOS  PALSAR  25-meter  Mosaic  Step-by-step  Manual  on  Extraction  of 
Forest  Cover  and  Change  Detection  Analysis  -  A  Manual  for  Users”.  Deutsche  Gesellschaft  fur 
Internationale  Zusammenarbeit  (GIZ)  GmbH. 

- .  (2013B):  “Support  to  Processing  of  Remote  Sensing  Data  for  the  Establishment  of  a  Pilot  MRV 

System  for  REDD+  on  Leyte  Island”.  GIZ  German  Technical  Cooperation  with  the  Philippines. 

Green,  A.A.,  M.  Berman,  P  Switzer  &  M.D.  Craig  (1988):  A  Transformation  for  Ordering  Multispectral 
Data  in  Terms  of  Image  Quality  with  Implications  for  Noise  Removal.  —IEEE  Transactions  on 
Geoscience  and  Remote  Sensing  26,  1,  65-74. 

GIZ  (2012):  Support  of  processing  Remote  Sensing  Data  for  the  Establishment  of  a  pilot  MRV  System  for 
REDD+  on  Leyte  Island.  Terms  of  Reference.  Gesellschaft  fur  Internationale  Zusammenarbeit  GIZ 
GmbH. 

Helmer,  E.H.  &  B.  Ruefenacht  (2007):  A  Comparison  of  radiometric  Normalization  Methods  when  filling 
Cloud  Gaps  in  Landsat  Imagery.  -Canadian  Journal  of  Remote  Sensing  33,  4,  325-340. 

Hsu,  C.-W.,  C.-C.  Chang  &  C.-J.  Lin  (2007):  A  practical  Guide  to  Support  Vector  Classification.  National 
Taiwan  University,  Department  of  Computer  Science  &  Information  Engineering  NTU  CSIE 
online.  <http:  / / ntu.  csie.ntu.edu.  tw /  -  cjlin /  papers/guide/ guide.pdf. 

Irish,  R.R.,  J.L.  Barker,  S.N.  Goward  &  T.  Arvidson  (2006):  Characterization  of  the  Landsat-7  ETM+ 
Automated  Cloud-Cover  Assessment  (ACCA)  Algorithm.  — Photogrammetric  Engineering  & 
Remote  Sensing  72,  10,  1179-1188. 

Kaufman,  Y.J.,  A.E.  Wald,  L.A.  Remer,  B.-C.  Gao,  R.-R.  Li  &  L.  Flynn  (1997):  The  MODIS  2.1  pm 
Channel  —  Correlation  with  Visible  Reflectance  for  Use  in  Remote  Sensing  of  Aerosol.  —IEEE 
Transactions  on  Geoscience  and  Remote  Sensing  35,  5,  1286-1298. 

Kistler,  R.,  E.  Kalnay,  W.  Collins,  S.  Saha,  G.  White,  J.  Woollen,  M.  Chelliah,  W.  Ebisuzaki,  M.  Kanamitsu, 
V.  Kousky,  H.  VanDeDool,  R.  Jenne  &  M.  Fiorino  (2001):  The  NCEP-NCAR  50-Year  Reanalysis: 
Monthly  Means  CD-ROM  and  Documentation.  —Bulletin  of  the  American  Meteorological  Society 
82,  2,  247-268. 


Kotchenova,  S.Y.,  E.F.  Vermote,  R.  Matarrese  &  F.  Klemm  (2006):  Validation  of  a  new  Vector  Version  of 
the  6S  radiative  Transfer  Code  for  atmospheric  Correction  ofMODIS  Data:  Part  I  -  Path  Radiance. 
-Applied  Optics  45,  26,  6762-6774. 

Malila  (1980):  Change  Vector  Analysis:  an  Approach  detecting  Forest  Changes  with  Fandsat.  -Proceedings 
of  the  6th  International  Symposium  on  Machine  Processing  of  Remotely  Sensed  Data,  326-335. 

Masek,  J.G.,  E.F  Vermote,  N.  Saleous,  R.  Wolfe,  F.G.  Hall,  F.  Huemmrich,  F.  Gao,  J.  Kutler,  andT.K.  Fim 
(2006):  A  Fandsat  Surface  Reflectance  Dataset  for  North  America,  1990-2000.  -IEEE  Geoscience 
and  Remote  Sensing  Fetters  3,  1,  68-72. 

Masek,  J.G.,  E.F  Vermote,  N.  Saleous,  R.  Wolfe,  F.G.  Hall,  F.  Huemmrich,  F.  Gao,  J.  Kutler,  andT.K.  Fim 
(2013):  FEDAPS  Calibration,  Reflectance,  Atmospheric  Correction  Preprocessing  Code,  Version  2. 
Model  product.  Available  online  from  Oak  Ridge  National  Faboratory  Distributed  Active  Archive 
Center,  Oak  Ridge,  Tennessee,  U.S.A.  <http://dx.doi.org/10.3334/  ORNFDAAC/1 146>. 

McPeters,  R.D.,  P.K.  Bhartia,  A.J.  Krueger,  J.R.  Herman,  C.G.  Wellemeyer,  C.J.  Seftor,  G.  Jaross,  O.  Torres, 
F.  Moy,  G.  Fabow,  W.  Byerly,  S.F.  Taylor,  T.  Swissler  &  R.P.  Cebula  (1998):  Earth  Probe  Total 
Ozone  Mapping  Spectrometer  (TOMS)  Data  Products  User’s  Guide.  NASA  Technical  Publication 
1998-206895.  Greenbelt,  Goddard  Space  Flight  Center. 

NASA  (2011):  Fandsat  7  Science  Data  Users  Handbook.  <http://landsathandbook.  gsfc.nasa.gov>. 

NCEP/NCAR  (2001):  NCEP/NCAR  Global  Reanalysis  Products,  1948-continuing.  National  Oceanic 
and  Atmospheric  Administration  (NOAA),  National  Weather  Service  (NWS),  National  Centers  for 
Environmental  Prediction  (NCEP).  —Research  Data  Archive  at  the  National  Center  for  Atmospheric 
Research  (NCAR),  Computational  and  Information  Systems  Faboratory  <http://rda.ucar.edu/ 
datasets/  ds090.0>. 

Nielsen,  A.A.  (2007):  The  Regularized  Iteratively  Reweighted  MAD  Method  for  Change  Detection  in 
Multi-  and  Hyperspectral  Data.  -IEEE  Transactions  on  Image  Processing  16,  2,  463-478. 

Richter,  R.,  T.  Kellenberger  &  H.  Kaufmann  (2009):  Comparison  of  Topographic  Correction  Methods. 
—Remote  Sensing  1,  184-196. 

Seifert-Granzin,  J.  (2013):  Conceptual  approach  to  REDD+  MRV  in  the  Philippines  -  An  Overview. 
Updated  Version  2.0.  Mesa  consult. 

Teillet,  P.M.,  B.  Guindon  &  D.G.  Goodenough  (1982):  On  the  Slope-Aspect  Correction  of  multispectral 
Scanner  Data.  -Canadian  Journal  of  Remote  Sensing  8,  84-106. 

Townshend,  J.R.,  J.G.  Masek,  C.  Huang,  E.F  Vermote,  F.  Gao,  S.  Channan,  J.O.  Sexton,  M.  Feng,  R. 
Narasimhan,  D.  Kim,  K.  Song,  D.  Song,  X.-P.  Song,  P.  Noojipady,  B.  Tan,  M.C.  Hansen,  M.  Fi  & 
R.E.  Wolfe  (2012):  Global  Characterization  of  Forest  Cover  using  Fandsat  Data:  Opportunities  and 
Challenges.  -International  Journal  of  Digital  Earth,  1-25. 

USGS  EROS  (2012):  GTOPO30.  Official  homepage.  United  States  Geological  Survey,  Earth  Resources 
Observation  and  Science.  <http://eros.usgs.gOv/#/Find_Data/  Products_and_Data_Available/ 
gtopo30_info>.  (Status:  2012-09-27)  (Access:  2013-06-06). 

Vermote,  E.F  &  N.  Saleous  (2007):  FEDAPS  Surface  Reflectance  Product  Description.  Version  2.  <http:// 
ledaps.nascom.nasa.gov>. 


Vermote,  E.  &  C.  Justice  (2010):  A  Surface  Reflectance  standard  Product  from  LDCM  and  supporting 
Activities.  June  2010  Update.  USGS  online  presence  <http://landsat.usgs.gov/documents/>  (status: 
28.12.2010)  (access:  14.08.2013). 

Vermote,  E.F.,  N.  El  Saleous,  C.O.  Justice,  Y.J.  Kaufman,  J.L.  Privette,  L.  Remer,  J.C.  Roger  &  D.  Tantre 
(1997):  Atmospheric  Correction  of  visible  to  middle-infrared  EOS-MODIS  data  over  land  surfaces: 
Background,  operational  algorithm  and  validation.  -Journal  of  Geophysical  Research  102,  17131- 
17141. 


Annexes 


Annexes 


Tab.  1:  List  of  integrated  images.  Separating  lines  delineate  the  four  produced  image  mosaics. 
One  full  island  coverage  per  acquisition  date  was,  if  possible,  produced  from  two  single  images. 
Master  acquisitions  for  radiometric  normalization  and  image  fusion  are  further  highlighted  in  bold 
letters. 
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Fig.  1:  Classification  result  of  1989-1990  image  mosaic. 
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Fig.  2:  Classification  result  of  1993-1994  image  mosaic. 
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Fig.  3:  Classification  result  of  1995-1996  image  mosaic. 
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Fig.  4:  Classification  result  of  2000-2001  image  mosaic. 
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